Testing Navier-Stokes with a moving body in package ViscousFlow

In [1]:

    

using ViscousFlow


    

    




┌ Info: Recompiling stale cache file /Users/jeff/.julia/compiled/v1.0/ViscousFlow/2kdbn.ji for ViscousFlow [103da179-b3e4-57c1-99a4-586354eb2c5a]
└ @ Base loading.jl:1190
[ Info: Building and caching LGF table

In [2]:

    

using Plots
pyplot()
clibrary(:colorbrewer)
default(grid = false)


    

    




┌ Warning: `getindex(o::PyObject, s::AbstractString)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o."s"` instead of `o["s"]`.
│   caller = top-level scope at none:0
└ @ Core none:0

In [3]:

    

Re = 200; # Reynolds number
U = 0.0; # Free stream velocity
U∞ = (U,0.0);


    

In [4]:

    

n = 51;
a = 0.5; b = 0.1;
body = Bodies.Plate(1.0,n)


    

    
Out[4]:





Plate with 51 points and length 1.0 and thickness 0.0
   Current position: (0.0,0.0)
   Current angle (rad): 0.0

In [6]:

    

plot(body)


    

    
Out[6]:

In [7]:

    

a = 0.25 # location of pitch axis, a = 0.5 is leading edge
ϕ = -π/2  # phase lag of pitch to heave
A = 0.25  # amplitude/chord
fstar = 1/π # fc/U
α₀ = 0 # mean angle of attack
Δα = 0.0 #10π/180 # amplitude of pitching
U₀ = 0.0 # translational motion (set to zero in place of free stream)

K = π*fstar # reduced frequency, K = πfc/U

oscil = RigidBodyMotions.PitchHeave(U₀,a,K,ϕ,α₀,Δα,A);
motion = RigidBodyMotion(oscil)


    

    
Out[7]:





Rigid Body Motion:
  ċ = 0.0 + 0.5im
  c̈ = 0.0 + 0.0im
  α̇ = 0.0
  α̈ = -0.0
  Oscillatory pitch-heave kinematics with
     Reduced frequency K = 1.0
     Heaving amplitude A = 0.25
     Pitching amplitude Δα = 0.0
     Pitch-to-heave lag ϕ = -1.5707963267948966

In [60]:

    

motion(0.0)


    

    
Out[60]:





(-0.25 + 0.0im, 0.0 + 0.5im, 0.0 + 0.0im, 0.0, 0.0, -0.0)

In [47]:

    

cent = 0.0 + 0.0im
α = 0.0
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration


    

    
Out[47]:





Plate with 51 points and length 1.0 and thickness 0.0
   Current position: (0.0,0.0)
   Current angle (rad): 0.0

In [49]:

    

plot(body)


    

    
Out[49]:

In [11]:

    

xlim = (-1.0,1.0)
ylim = (-1.0,1.0)


    

    
Out[11]:





(-1.0, 1.0)

In [50]:

    

X = VectorData(body.x,body.y);
X̃ = VectorData(body.x̃,body.ỹ);


    

In [53]:

    

Δx = 0.02;
Δt = min(0.5*Δx,0.5*Δx^2*Re);


    

In [63]:

    

sys = NavierStokes(Re,Δx,xlim,ylim,Δt,U∞ = U∞, X̃ = X̃, isstore = true, isstatic = false)


    

    
Out[63]:





Navier-Stokes system on a grid of size 102 x 102

In [64]:

    

w₀ = Nodes(Dual,size(sys));
u = (w₀,[real(cent),imag(cent),α])
xg, yg = coordinates(w₀,dx=Δx,I0=Systems.origin(sys))
f = (VectorData(X̃),Vector{Float64}());


    

In [65]:

    

using LinearAlgebra
plans = ((t,u) -> Fields.plan_intfact(t,u,sys),(t,u)-> I)


    

    
Out[65]:





(getfield(Main, Symbol("##81#83"))(), getfield(Main, Symbol("##82#84"))())

In [80]:

    

function TimeMarching.plan_constraints(u::Tuple{Nodes{Dual,NX,NY},Vector{Float64}},t,sys::NavierStokes{NX,NY,N,false}) where {NX,NY,N}

  # for now, just assume that there is only one body. will fix later.

  xc, yc, α = u[2]
  T = Bodies.RigidTransform((xc,yc),α)
  # should be able to save some time and memory allocation here...
  x, y = T(sys.X̃.u,sys.X̃.v)
  X = VectorData(x,y)
  regop = Regularize(X,cellsize(sys);issymmetric=true,I0=origin(sys))
  if sys._isstore
    Hmat, Emat = RegularizationMatrix(regop,VectorData{N}(),Edges{Primal,NX,NY}())
      
    sys.Hmat = Hmat
    sys.Emat = Emat
    
    return (f->TimeMarching.B₁ᵀ(f,sys), f->zeros(Float64,size(u[2]))),
           (w->TimeMarching.B₂(w,sys), u->Vector{Float64}())
  else
    return (f->TimeMarching.B₁ᵀ(f,regop,sys), f->zeros(Float64,size(u[2]))),
           (w->TimeMarching.B₂(w,regop,sys), u->Vector{Float64}())
  end


end


    

In [82]:

    

ifherk = IFHERK(u,f,sys.Δt,
                plans,
                (u,t) -> TimeMarching.plan_constraints(u,t,sys),
                ((u,t) -> Systems.r₁(u,t,sys,motion),
                 (u,t) -> Systems.r₂(u,t,sys,motion)),tol=1e-3,rk=TimeMarching.RK31,isstored=true,isstaticconstraints=false)


    

    
Out[82]:





Order-3 IF-HERK integrator with
   State of type Tuple{Nodes{Dual,102,102},Array{Float64,1}}
   Force of type Tuple{VectorData{51},Array{Float64,1}}
   Time step size 0.01

In [83]:

    

t = 0.0
w₀ .= 0.0
u = (deepcopy(w₀),[real(cent),imag(cent),α])
f = (VectorData(X̃),Vector{Float64}());
tf = 0.1;
T = Δt:Δt:tf;
fx = Float64[];
fy = Float64[];
thist = [];
uhist = [];
tsample = 0.1;


    

In [86]:

    

@time for ti in T
    global t, u, f = ifherk(t,u)
    
    # save data for later use
    push!(thist,t)
    push!(fx,sum(f[1].u)*Δx^2)
    push!(fy,sum(f[1].v)*Δx^2)
    (isapprox(mod(t,tsample),0,atol=1e-12) || isapprox(mod(t,tsample),tsample,atol=1e-12)) ? push!(uhist,u) : nothing
end
println("solution completed through time t = ",t)


    

    




 19.901916 seconds (189.39 M allocations: 6.132 GiB, 7.45% gc time)
solution completed through time t = 0.20000000000000004

In [87]:

    

plot(xg,yg,uhist[1][1],levels=range(-1,1,length=30),color = :RdBu)
Tr = RigidTransform(uhist[1][2])
Tr(body) # transform the body to the current configuration
plot!(body,linecolor=:black)


    

    
Out[87]:

In [88]:

    

px = plot(thist,2*fx,ylim=(-3,3),xlabel="Convective time",ylabel="\$C_D\$",legend=false)
py = plot(thist,2*fy,ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=false)
plot(px,py)


    

    
Out[88]:

In [90]:

    

ψ = deepcopy(u[1])
ψ .= -(sys.L\u[1])*Δx .+ repeat(collect(yg)',size(sys,1),1)
plot(xg,yg,ψ,trim=1,clim=(-0.1,0.1))
plot!(body,fillcolor=:black,linecolor=:black)


    

    
Out[90]:

In [ ]: